Craugastor Ranoides Microbiome Analysis
Xav Harrison
2025-07-29
- 0.1 Libraries
- 0.2 Global Plotting Options
- 0.3 Our Data
- 0.4 Frog Metadata
- 0.5 FILTER UNWANTED TAXONOMY
- 0.6 Identifying and Removing Decontaminants
- 0.7 Filter Based on Prevalence and Abundance
- 0.8 Mock Community
- 0.9 Assessing Post-QC Sample Coverage and Library Sizes
- 0.10 Reads by Site and Study
- 1 ANALYSIS
- 2 GLLVM
- 3 UpSet Plot
- 4 Sample Similarity Network
- 5 TEMPORAL ANALYSIS
- 6 Temporal GLLVM
- 7 FUNCTIONAL ANALYSIS
0.1 Libraries
Packages we’ll need for all the stats
library(phyloseq) #microbiome data handling
library(ggplot2) #plotting
library(dplyr) # data handling ##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union library(vegan) #plotting of community data like bacterial microbiomes and statistical models## Loading required package: permute library(decontam) #identifying contaminants in our sequencing workflow
library(microbiome) # microbiome plotting and convenience functions##
## microbiome R package (microbiome.github.com)
##
##
##
## Copyright (C) 2011-2022 Leo Lahti,
## Sudarshan Shetty et al. <microbiome.github.io>##
## Attaching package: 'microbiome'## The following object is masked from 'package:vegan':
##
## diversity## The following object is masked from 'package:ggplot2':
##
## alpha## The following object is masked from 'package:base':
##
## transform library(ggordiplots) # tidyverse versions of the vegan plots ## Loading required package: glue library(MASS) #generalized linear models of pathogen infection data (Negative Binomial Error Structure)##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select library(pheatmap) #pretty heatmaps
library(RColorBrewer) #Add Custom Colours (e.g. Colour Blind Friendy)
library(MuMIn)
#devtools::install_github("BlakeRMills/MoMAColors")
library(MoMAColors)
library(tidybayes)
library(cowplot)
library(ggvenn) #install.packages('ggvenn')## Loading required package: grid library(janitor)##
## Attaching package: 'janitor'## The following objects are masked from 'package:stats':
##
## chisq.test, fisher.test library(brms)## Loading required package: Rcpp## Loading 'brms' package (version 2.22.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').##
## Attaching package: 'brms'## The following objects are masked from 'package:tidybayes':
##
## dstudent_t, pstudent_t, qstudent_t, rstudent_t## The following object is masked from 'package:MuMIn':
##
## loo## The following object is masked from 'package:phyloseq':
##
## nsamples## The following object is masked from 'package:stats':
##
## ar library(microbiomeutilities)## Warning: replacing previous import 'ggplot2::alpha' by 'microbiome::alpha' when
## loading 'microbiomeutilities'## Registered S3 method overwritten by 'broom':
## method from
## nobs.multinom MuMIn##
## Attaching package: 'microbiomeutilities'## The following object is masked from 'package:microbiome':
##
## add_refseq library(gllvm)## Loading required package: TMB##
## Attaching package: 'gllvm'## The following objects are masked from 'package:MuMIn':
##
## AICc, coefplot## The following object is masked from 'package:vegan':
##
## ordiplot## The following object is masked from 'package:stats':
##
## simulate library(ComplexUpset)0.2 Global Plotting Options
Quick access options for output graphics to make legends and axis labels larger / more legible
#Global Plot Options
plotopts<- theme(axis.text=element_text(size=20),axis.title=element_text(size=22),strip.text=element_text(size=22),legend.title = element_text(size=20),legend.text = element_text(size=20))
## Global Site Colors
#sitecols<-c(brewer.pal(8,"Set2"),brewer.pal(9,"Paired")[9])
###### Stop Warnings turning into Errors
# options(warn=1)0.3 Our Data
Load in Our Phyloseq Object
load('CostaRica_Ranoides_Phyloseq.RData')0.4 Frog Metadata
frogmeta<-read.csv('Ranoides_metadata_newsites_oct24.csv',header=T)
head(frogmeta)## Plate Well Sample Library Protocol PCR Barcode.Name Barcode.Sequence
## 1 1 A1 RP01 RP01-l1 Other none UDP0193 TCTCATGATA-AACACGTGGA
## 2 1 B1 RP02 RP02-l1 Other none UDP0194 CGAGGCCAAG-GTGTTACCGG
## 3 1 C1 RP03 RP03-l1 Other none UDP0195 TTCACGAGAC-AGATTGTTAC
## 4 1 D1 RP04 RP04-l1 Other none UDP0196 GCGTGGATGG-TTGACCAATG
## 5 1 E1 RP05 RP05-l1 Other none UDP0197 TCCTGGTTGT-CTGACCGGCA
## 6 1 F1 RP06 RP06-l1 Other none UDP0198 TAATTCTGCT-TCTCATCAAT
## Conc Flags class SV..mm. Site.name GPS.code Lat Long
## 1 NA NA 41.0 Brazo del potero grande CI01 10.88052 -85.69252
## 2 NA NA 35.0 Brazo del potero grande CI01 10.88052 -85.69252
## 3 NA NA 26.2 Brazo del potero grande CI01 10.88052 -85.69252
## 4 NA NA 20.0 Brazo del potero grande CI02 10.89880 -85.75618
## 5 NA NA 28.3 Brazo del potero grande CI03 10.87838 -85.69412
## 6 NA NA female 39.0 Brazo del potero grande CI03 10.87838 -85.69412
## M.A.S.L...GPS. Date Time Swabbed Recaptur..New KH.note pH_mean pH_sd
## 1 372.6681 04/01/2022 19:05 Yes N 7.87 0.29
## 2 372.6681 04/01/2022 19:05 Yes N 7.87 0.29
## 3 372.6681 04/01/2022 19:05 Yes N 7.87 0.29
## 4 85.4131 04/01/2022 20:01 Yes N 7.87 0.29
## 5 333.6929 04/01/2022 20:32 Yes N 7.87 0.29
## 6 333.6929 04/01/2022 20:32 Yes N 7.87 0.29
## pH_n pH.FP_mean pH.FP_sd pH.FP_n Site
## 1 12 8.1 0.19 4 BdPG-B
## 2 12 8.1 0.19 4 BdPG-B
## 3 12 8.1 0.19 4 BdPG-B
## 4 12 8.1 0.19 4 BdPG-B
## 5 12 8.1 0.19 4 BdPG-A
## 6 12 8.1 0.19 4 BdPG-A
## Withinsite.options
## 1 Split.site_Between.upper.and.lower.waterfalls
## 2 Split.site_Between.upper.and.lower.waterfalls
## 3 Split.site_Between.upper.and.lower.waterfalls
## 4 Split.site_Between.upper.and.lower.waterfalls
## 5 Split.site_After.lower.waterfall
## 6 Split.site_After.lower.waterfall
## Note.on.withinsite.split WS.pH_mean WS.pH_sd WS.pH_n WS.pH.FP_mean
## 1 based.on.field.note.site.ecology 7.88 0.21 6 7.98
## 2 based.on.field.note.site.ecology 7.88 0.21 6 7.98
## 3 based.on.field.note.site.ecology 7.88 0.21 6 7.98
## 4 based.on.field.note.site.ecology 7.88 0.21 6 7.98
## 5 based.on.field.note.site.ecology 7.98 0.27 5 8.14
## 6 based.on.field.note.site.ecology 7.98 0.27 5 8.14
## WS.pH.FP_sd WS.pH.FP_n PIT.tag.code Elastomer. Colour.left Colour.right
## 1 NA 1 yes orange orange
## 2 NA 1 yes red red
## 3 NA 1 yes red orange
## 4 NA 1 yes
## 5 0.21 3 yes orange red
## 6 0.21 3 3D6.15340D80C8 No
## Notes
## 1 From Cerro el Ingles
## 2
## 3
## 4 too small to mark
## 5
## 6
## People.on.trip
## 1 Harrison, Herborn, North, Wilson, Puschendorf, Siegfried
## 2 Harrison, Herborn, North, Wilson, Puschendorf, Siegfried
## 3 Harrison, Herborn, North, Wilson, Puschendorf, Siegfried
## 4 Harrison, Herborn, North, Wilson, Puschendorf, Siegfried
## 5 Harrison, Herborn, North, Wilson, Puschendorf, Siegfried
## 6 Harrison, Herborn, North, Wilson, Puschendorf, Siegfriedfrogmeta<-clean_names(frogmeta)
rownames(frogmeta)<-frogmeta$sample
########### CLEAN UP NEGATIVES AND SAMPLE TYPE VARIABLES
frogmeta$sampletype<-"frog"
frogmeta$sampletype[grep("NEG|NTC",frogmeta$sample)]<-"Negative_Control"
frogmeta$sampletype[grep("MOCK|Mock",frogmeta$sample)]<-"Positive_Control"
########### DATE MANIPULATION
table(frogmeta$date_time)## < table of extent 0 >########## CODE UP TRIPS
frogmeta$study<-"spatial"
frogmeta$study[which(frogmeta$sampletype %in% c("Negative_Control","Positive_Control"))]<-"controls"
frogmeta$study[which(frogmeta$site_name=="" & !frogmeta$sampletype %in% c("Negative_Control","Positive_Control"))]<-"temporal"
table(frogmeta$study)##
## controls spatial temporal
## 6 48 23##Code Up Timepoints for Temporal
frogmeta$timepoint<-NA
frogmeta$timepoint[frogmeta$study=="temporal"]<-"PreTranslocation"
frogmeta$timepoint[frogmeta$study=="temporal" & grepl("_",frogmeta$sample)]<-"PostTranslocation"
##Code Up Indivieduals for Temporal
frogmeta$id<-NA
id_temp<- as.numeric(gsub(".*?([0-9]+).*", "\\1", frogmeta$sample)) ## Warning: NAs introduced by coercion frogmeta$id[frogmeta$study=="temporal"]<-id_temp[frogmeta$study=="temporal"]
table(frogmeta$id)##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14
## 1 2 2 2 1 2 2 1 2 2 1 1 2 2############# STITCH TOGETHER
sample_data(ps)<-frogmeta
ps## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 31120 taxa and 76 samples ]
## sample_data() Sample Data: [ 76 samples by 47 sample variables ]
## tax_table() Taxonomy Table: [ 31120 taxa by 6 taxonomic ranks ]
## refseq() DNAStringSet: [ 31120 reference sequences ] #write.csv(frogmeta,'Ranoides Sample Metadata Final.csv')0.5 FILTER UNWANTED TAXONOMY
#Prune Taxa With No Phylum Assignment
ps_prune<-prune_taxa(as.vector(!is.na(tax_table(ps)[,2])),ps)
ntaxa(ps)-ntaxa(ps_prune) #326## [1] 326 #Prune Chloroplasts
ps_prune_nochloro<-prune_taxa(as.vector(tax_table(ps_prune)[,4]!="Chloroplast"),ps_prune)
ntaxa(ps_prune)-ntaxa(ps_prune_nochloro) #2827## [1] 2827 #Remove Mitochondria
ps_prune_nochloro_nomito<-prune_taxa(as.vector(tax_table(ps_prune_nochloro)[,5]!="Mitochondria"),ps_prune_nochloro)
#Remove
ps_prune_nochloro_nomito_noarchaea<-prune_taxa(as.vector(tax_table(ps_prune_nochloro_nomito)[,1]!="Archea"),ps_prune_nochloro_nomito)
ps_prune_nochloro_nomito_noarchaea## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 23736 taxa and 76 samples ]
## sample_data() Sample Data: [ 76 samples by 47 sample variables ]
## tax_table() Taxonomy Table: [ 23736 taxa by 6 taxonomic ranks ]
## refseq() DNAStringSet: [ 23736 reference sequences ]0.6 Identifying and Removing Decontaminants
## Flag Negatives
sample_data(ps_prune_nochloro_nomito_noarchaea)$is_negative<- sample_data(ps_prune_nochloro_nomito_noarchaea)$sampletype=="Negative_Control"
##Inspect Library SIzes
df <- as.data.frame(sample_data(ps_prune_nochloro_nomito_noarchaea)) # Put sample_data into a ggplot-friendly data.frame
df$LibrarySize <- sample_sums(ps_prune_nochloro_nomito_noarchaea)
df <- df[order(df$LibrarySize),]
df$Index <- seq(nrow(df))
ggplot(data=df, aes(x=Index, y=LibrarySize, color=is_negative)) + geom_point()
#Negatives Based on Prevalence // Threshold 0.5
contamdf.prev <- isContaminant(ps_prune_nochloro_nomito_noarchaea, method="prevalence", neg="is_negative",threshold=0.5)
table(contamdf.prev$contaminant) #51##
## FALSE TRUE
## 23685 51#Plot
ps.pa <- transform_sample_counts(ps_prune_nochloro_nomito_noarchaea, function(abund) 1*(abund>0))
ps.pa.neg <- prune_samples(sample_data(ps.pa)$is_negative == TRUE, ps.pa)
ps.pa.pos <- prune_samples(sample_data(ps.pa)$is_negative == FALSE, ps.pa)
# Make data.frame of prevalence in positive and negative samples
df.pa <- data.frame(pa.pos=taxa_sums(ps.pa.pos), pa.neg=taxa_sums(ps.pa.neg),
contaminant=contamdf.prev$contaminant)
ggplot(data=df.pa, aes(x=pa.neg, y=pa.pos, color=contaminant)) + geom_point() +
xlab("Prevalence (Negative Controls)") + ylab("Prevalence (True Samples)")
####### PRUNE OUT
ps_clean<-prune_taxa(contamdf.prev$contaminant==FALSE,ps_prune_nochloro_nomito_noarchaea)
######## Inspect POSITIVES
ps_clean_positives<-prune_samples(sample_data(ps_clean)$sampletype=="Positive_Control",ps_clean)
ps_clean_positives<-filter_taxa(ps_clean_positives, function(x) (sum(x) > 0),TRUE)
data.frame(reads=taxa_sums(ps_clean_positives),genus=tax_table(ps_clean_positives)[,6])## reads Genus
## a4cbe987942964b584e95e4efab6a176 246898 Bacillus
## 8ae518dbb29595b3f79214be0b589066 39343 Listeria
## d46e2205f0c6ecf67b51f83d111c509c 127844 Escherichia-Shigella
## d829bee4984f82ffc2453212157caf96 5 Bradyrhizobium
## 4cbfff144d4e7a4e0f4619ed505be070 97729 Salmonella
## ff9d93d7b7e46787568f2d241caeaf3b 54763 Pseudomonas
## 65d43491988bfe557da4d86a5ba25dae 20271 Staphylococcus
## b92c7cfb137fee89c86d0a74209d501a 44 <NA>
## 3e86e9620a4aceb493a34d5c47723513 15866 Salmonella
## 4dbd335b510606def06a34bef91c5214 20 Acetobacter
## db29dfd5db5e2501ed9deadabc7dd91d 23 Acetobacter
## 892a0fc53b26d6a2c193d145b2464da7 9 Chryseobacterium
## 86748e6e46eddb955f6918863b3cf191 10 <NA>
## bdf8a26094624622d68509a87fa75ba7 21 Bacillus
## ac3e504fcfb796da3a3c303eb96b4999 4797 Escherichia-Shigella
## 674d3fca8b60d684ea44b52ad3fd4c73 9 Polaromonas
## 8df399e29acb5a583f449a82db4e17c9 13 Lactiplantibacillus
## bd569194c3d7c9dc1817718bf17a19b2 13 Ligilactobacillus
## ce06233c4e2a93fd65d6a65637073148 1536 Bacillus
## be5cc85bd2ac16b6ae446fd24b1e2308 8 Turicibacter
## 569653b1659271a290facfcedd0de061 214 Enterococcus
## e266243eb28b5f38b6f567e2603d27f3 5 Bifidobacterium
## 1ed490b0a488876bb6f2f1167db28782 23 <NA>
## 7b08fffde750ca1296315e18f17ff004 271 <NA>
## 1309d2aeee2812f04dd6acf447d53af9 88 Bacillus
## 8e850704e6d34fe2ba5cd32b2f184224 108 Listeria
## c27583bfe4cd656b824bf5d33c8a9827 8 <NA>
## a3f9df4186a4f00562674fb08daca5d4 89 Paenibacillus
## fd44d4cb468fd7dc9b3227867714ed87 4 Bacteroides
## 48322a34013aed7ee00b201106330273 3 Faecalibacterium
## b9fcd7d71b74853248517b892d03a745 8 Flavobacterium0.7 Filter Based on Prevalence and Abundance
####### Read Threshold and Sample Prevalence Threshold
ps_clean_filter = filter_taxa(ps_clean, function(x) sum(x > 20) > (5), TRUE)
######### SUBSET TO ONLY FROGS AND WATER (ALSO STRIP NAs)
ps_clean_samples<-prune_samples(sample_data(ps_clean_filter)$sampletype %in% c("frog"),ps_clean_filter)
ps_clean_samples## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 1038 taxa and 70 samples ]
## sample_data() Sample Data: [ 70 samples by 48 sample variables ]
## tax_table() Taxonomy Table: [ 1038 taxa by 6 taxonomic ranks ]
## refseq() DNAStringSet: [ 1038 reference sequences ]0.8 Mock Community
##Inspect Mock After Removing ASVs filtered at previous step
ps_mock<-prune_samples(sample_data(ps_clean_filter)$sampletype %in% c("Positive_Control"),ps_clean_filter)
ps_mock<-filter_taxa(ps_mock, function(x) (sum(x) > 0),TRUE)
mock_dataframe<-data.frame(reads=taxa_sums(ps_mock),genus=tax_table(ps_mock)[,6])
#Imposters in Mocks
mock_imposters<-rownames(mock_dataframe)[mock_dataframe[,1]<1000]
#What Proportion of Dataset
taxa_sums_data<-data.frame(taxa_sums(ps_clean_filter))
subset(taxa_sums_data,rownames(taxa_sums_data) %in% mock_imposters)## taxa_sums.ps_clean_filter.
## d829bee4984f82ffc2453212157caf96 137873
## b92c7cfb137fee89c86d0a74209d501a 43082
## 4dbd335b510606def06a34bef91c5214 14180
## db29dfd5db5e2501ed9deadabc7dd91d 13502
## 892a0fc53b26d6a2c193d145b2464da7 9775
## 86748e6e46eddb955f6918863b3cf191 9693
## bdf8a26094624622d68509a87fa75ba7 7046
## 674d3fca8b60d684ea44b52ad3fd4c73 4701
## 8df399e29acb5a583f449a82db4e17c9 4207
## bd569194c3d7c9dc1817718bf17a19b2 3500
## be5cc85bd2ac16b6ae446fd24b1e2308 1407
## 569653b1659271a290facfcedd0de061 982
## e266243eb28b5f38b6f567e2603d27f3 9510.9 Assessing Post-QC Sample Coverage and Library Sizes
Now that we’ve done that, we can check what our post-QC library sizes are. It’s a good idea to report these in manuscripts and thesis chapters. Here was ask what the minimum, maximum and mean library sizes-per-samples are.
############## POST QC LIBRARY SIZES
###############
# Post-QC Library Stats
###############
mean(sample_sums(ps_clean_samples))## [1] 103222.5 range(sample_sums(ps_clean_samples))## [1] 15186 1755580.10 Reads by Site and Study
#Make a data frame of read depths
frog_reads<-data.frame(reads=sample_sums(ps_clean_samples))
range(frog_reads$reads)## [1] 15186 175558 #Add on the sample ID
frog_reads$sample<-rownames(frog_reads)
#Join on the Metadata
frog_meta<-as(sample_data(ps_clean_samples),"data.frame")
frog_reads<-left_join(frog_reads,frog_meta,"sample")
#Some Boxplots of Coverage by NEW Population using ggplot2
ggplot(frog_reads,aes(x=site,y=reads)) + geom_violin()
#Some Boxplots of Coverage by OLD Population using ggplot2
ggplot(frog_reads,aes(x=site_name,y=reads)) + geom_violin()
#Some Boxplots of Coverage by Population using ggplot2
ggplot(frog_reads,aes(x=study,y=reads)) + geom_violin()
1 ANALYSIS
1.1 Rarefaction Curves
#Strip Out the OTU Table
ps_clean_otutable<-as(t(otu_table(ps_clean_samples)),"matrix")
class(ps_clean_otutable)<-"matrix"
#Rarefaction Curves
raremax <- min(rowSums(ps_clean_otutable))
ps_clean_otutable[1,1]<-1
rarecurve(ps_clean_otutable,step=50,sample = 100000)
1.2 Rarefying Samples
#remind ourselves what the minimum coverage is
min(sample_sums(ps_clean_samples))## [1] 15186 #Histogram
hist(sample_sums(ps_clean_samples),breaks=10000)
# Rarefy to lowest library size, set random number see for reproducibility
ps_rare<-rarefy_even_depth(ps_clean_samples,rngseed = 150517,sample.size=15000)## `set.seed(150517)` was used to initialize repeatable random subsampling.## Please record this for your records so others can reproduce.## Try `set.seed(150517); .Random.seed` for the full vector## ...1.2.1 Estimate Richness
#Estimate Observed RIchness and Shannon Diversity
frog_rich<-estimate_richness(ps_rare,measures=c('Shannon','Observed'))
frog_rich$sample<-gsub("^X","",rownames(frog_rich))
#Add On Sample metadata
frog_rich<-left_join(frog_rich,frog_meta,"sample")
## ID Replication?
table(paste0(frog_rich$colour_left,frog_rich$colour_right))##
## blueblue bluered blueyellow greengreen greenorange
## 50 1 1 1 1 1
## greenyellow orangeblue orangegreen orangeorange orangered orangeyellow
## 1 1 1 1 1 1
## red blue redorange redred redyellow yellowblue yellowgreen
## 1 1 1 1 1 1
## yelloworange yellowred yellowyellow
## 1 1 1## Site Colours
sitecols<-brewer.pal(7,"Paired")[c(1,3,5)]1.2.2 Does Richness Vary By Site - Spatial
#Spatial Project Only
frog_spatial<-subset(frog_rich,study=="spatial")
table(frog_spatial$site)##
## BdPG-A BdPG-B M-A M-B P-A P-B
## 7 6 3 8 6 18 frog_spatial$site_name2<-frog_spatial$site_name
frog_spatial$site_name_short<-sapply(strsplit(frog_spatial$site,"-"),"[",1)
frog_spatial$site_name2[frog_spatial$site_name2=="Brazo del potero grande"]<-"Brazo del \n potrero grande "
#Plot by Site OLD
site_richness1<- frog_spatial %>% ggplot(.,aes(y=Observed,x=site_name2)) + geom_point(aes(fill=site_name2),shape=21,size=6,colour="white") + stat_pointinterval(shape=21,point_fill="white",point_size=5,colour="black",slab_alpha=0.5) + theme_bw() + plotopts
site_richness2 <- site_richness1 + guides(fill="none") + labs(y="Alpha Diversity (Richness)",x="Site") + scale_fill_manual(values=sitecols) + theme(axis.text.x=element_text(angle=25,hjust=1))
site_richness2
#ggsave2('Ranoides Richness by Site.pdf',site_richness2,width=20,height=15,units="cm")
#Plot by Site NEW
site_new_richness1<- frog_spatial %>% ggplot(.,aes(y=Observed,x=site)) + geom_point(aes(fill=site),shape=21,size=6,colour="white",position = position_jitter(width=0.15)) + stat_pointinterval(shape=21,point_fill="white",point_size=5,colour="black",slab_alpha=0.5) + theme_bw() + plotopts
site_new_richness2 <- site_new_richness1 + guides(fill="none") + labs(y="Alpha Diversity \n(Richness)",x="Site") + scale_fill_brewer(palette = "Paired") + theme(axis.text.x=element_text(angle=25,hjust=1))
site_new_richness2
#Shannon by Site NEW
site_new_shannon1<- frog_spatial %>% ggplot(.,aes(y=Shannon,x=site)) + geom_point(aes(fill=site),shape=21,size=6,colour="white",position = position_jitter(width=0.15)) + stat_pointinterval(shape=21,point_fill="white",point_size=5,colour="black",slab_alpha=0.5) + theme_bw() + plotopts
site_new_shannon2 <- site_new_shannon1 + guides(fill="none") + labs(y="Alpha Diversity \n(Richness)",x="Site") + scale_fill_brewer(palette = "Paired") + theme(axis.text.x=element_text(angle=25,hjust=1))
site_new_shannon2 
#Plot by Size
svl_richness1<- frog_spatial %>% ggplot(.,aes(y=Observed,x=sv_mm)) + geom_smooth(method = "lm") + geom_point(aes(fill=site_name_short),shape=21,size=6,colour="white") + theme_bw() + plotopts
svl_richness2 <- svl_richness1 +theme(legend.position = "top") + labs(fill="Site",y="Alpha Diversity \n(Richness)",x="Snout Vent Length (mm)") + scale_fill_manual(values=sitecols) + facet_wrap(.~site_name_short)
svl_richness2## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
ggsave2('Ranoides Richness by SVL.pdf',svl_richness2,width=20,height=15,units="cm")## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 1 row containing non-finite outside the scale range (`stat_smooth()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).#RIchness by pH
ph_richness1<- frog_spatial %>% ggplot(.,aes(y=Observed,x=p_h_mean)) + geom_smooth(method = "lm") + geom_point(aes(fill=site_name_short),shape=21,size=6,colour="white") + theme_bw() + plotopts
ph_richness2 <- ph_richness1 +theme(legend.position = "top") + labs(fill="Site",y="Alpha Diversity \n(Richness)",x="pH") + scale_fill_manual(values=sitecols) #+ facet_wrap(.~site)
ph_richness2## `geom_smooth()` using formula = 'y ~ x'
ggsave2('pH Richness.png',ph_richness2,dpi=150)## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x' with(frog_spatial,cor.test(Observed,p_h_mean))##
## Pearson's product-moment correlation
##
## data: Observed and p_h_mean
## t = -1.0378, df = 46, p-value = 0.3048
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.4174496 0.1388462
## sample estimates:
## cor
## -0.15125571.3 Richness Models
richmod1<-brm(Observed ~ sv_mm + site_name_short,data=frog_spatial,family=negbinomial())## Warning: Rows containing NAs were excluded from the model.## Compiling Stan program...## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.5)’
## using SDK: ‘’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## 679 | #include <cmath>
## | ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1## Start sampling##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 2.3e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.23 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.037 seconds (Warm-up)
## Chain 1: 0.021 seconds (Sampling)
## Chain 1: 0.058 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 3e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.041 seconds (Warm-up)
## Chain 2: 0.022 seconds (Sampling)
## Chain 2: 0.063 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 3e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.038 seconds (Warm-up)
## Chain 3: 0.022 seconds (Sampling)
## Chain 3: 0.06 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 3e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.03 seconds (Warm-up)
## Chain 4: 0.024 seconds (Sampling)
## Chain 4: 0.054 seconds (Total)
## Chain 4:summary(richmod1)## Family: negbinomial
## Links: mu = log; shape = identity
## Formula: Observed ~ sv_mm + site_name_short
## Data: frog_spatial (Number of observations: 47)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 5.06 0.22 4.62 5.51 1.00 4276 2977
## sv_mm 0.01 0.01 -0.00 0.02 1.00 5192 3094
## site_name_shortM -0.05 0.17 -0.39 0.29 1.00 3530 2871
## site_name_shortP 0.16 0.14 -0.13 0.43 1.00 3541 2823
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 6.33 1.39 3.97 9.39 1.00 3711 2559
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).conditional_effects(richmod1)
#Add Null Model
frog_spatial_nona<-subset(frog_spatial,!is.na(sv_mm)) #remove NA value for SVL
richmod1_null<-brm(Observed ~ 1,data=frog_spatial_nona,family=negbinomial())## Compiling Stan program...
## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.5)’
## using SDK: ‘’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## 679 | #include <cmath>
## | ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1## Start sampling##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 1.8e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.18 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.024 seconds (Warm-up)
## Chain 1: 0.021 seconds (Sampling)
## Chain 1: 0.045 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 4e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.023 seconds (Warm-up)
## Chain 2: 0.021 seconds (Sampling)
## Chain 2: 0.044 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 4e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.022 seconds (Warm-up)
## Chain 3: 0.024 seconds (Sampling)
## Chain 3: 0.046 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 6e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.025 seconds (Warm-up)
## Chain 4: 0.025 seconds (Sampling)
## Chain 4: 0.05 seconds (Total)
## Chain 4:summary(richmod1_null)## Family: negbinomial
## Links: mu = log; shape = identity
## Formula: Observed ~ 1
## Data: frog_spatial_nona (Number of observations: 47)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 5.46 0.06 5.34 5.58 1.00 2955 2322
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 6.07 1.22 3.92 8.69 1.00 3076 2691
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).richmod1_null<-add_criterion(richmod1_null,"loo")
richmod1<-add_criterion(richmod1,"loo")
print(loo_compare(richmod1,richmod1_null),simplify = F)## elpd_diff se_diff elpd_loo se_elpd_loo p_loo se_p_loo looic
## richmod1_null 0.0 0.0 -280.0 3.8 1.6 0.3 560.0
## richmod1 -0.6 1.9 -280.6 4.1 4.4 0.8 561.2
## se_looic
## richmod1_null 7.6
## richmod1 8.1#Pull Out Predicted Data
rich_post<-conditional_effects(richmod1)[2]$site_name_short
### Posterior Plots
site_new_richness_bayes1<- ggplot() + geom_point(data=frog_spatial,aes(fill=site_name_short,y=Observed,x=site_name_short),shape=21,size=6,colour="white",position = position_jitter(width=0.15)) + geom_errorbar(data=rich_post,aes(x=site_name_short,ymin=lower__,ymax=upper__),width=0.15) + geom_point(data=rich_post,aes(x=site_name_short,y=estimate__),shape=21,size=5,fill="white") + theme_bw() + plotopts
site_new_richness_bayes2 <- site_new_richness_bayes1 + guides(fill="none") + labs(y="Alpha Diversity \n(Richness)",x="Site") + scale_fill_manual(values=sitecols) + theme(axis.text.x=element_text(angle=25,hjust=1))
site_new_richness_bayes2
1.4 ORDINATIONS
1.4.1 Centred-Log Ratio Ordinations (CLR)
Will Just do on Frogs for this now
### Just Frogs
ps_clean_spatial<-prune_samples(sample_data(ps_clean_samples)$study=="spatial",ps_clean_samples)
sample_data(ps_clean_spatial)$site_name_short<-sapply(strsplit(sample_data(ps_clean_spatial)$site,"-"),"[",1)
#CLR TRANSFORM
ps_spatial_clr<-microbiome::transform(ps_clean_spatial, 'clr')
#Ordinate
ord_spatial_clr <- phyloseq::ordinate(ps_spatial_clr, "RDA")
phyloseq::plot_ordination(ps_spatial_clr, ord_spatial_clr, type="samples", color="site_name")
#GG Version NEW SITES
gg_ordiplot(ord_spatial_clr,sample_data(ps_spatial_clr)$site,spiders = TRUE,ellipse = T)
#GG Version
gg_ordiplot(ord_spatial_clr,sample_data(ps_spatial_clr)$site_name_short,spiders = TRUE,ellipse = T)
#### Extract for Plotting
#Plot and Extract for GGPLot2
x1 <- gg_ordiplot(ord_spatial_clr,sample_data(ps_spatial_clr)$site_name_short,spiders = TRUE,ellipse = T,plot=F)
x1_plotdata<- x1$df_ord
x1_plotdata$Group<-as.character(x1_plotdata$Group)
#x1_plotdata$Group[which(x1_plotdata$Group=="Brazo del potero grande")]<-"Cerro Ingles"
### Site Labels
#Plot Spiders
frog_site_clr_plot1<-ggplot(x1_plotdata,aes(x=x,y=y)) + geom_segment(data=x1$df_spiders, aes(x=cntr.x, xend=x, y=cntr.y, yend=y, color=Group), show.legend = FALSE)
#Plot Points
frog_site_clr_plot2<- frog_site_clr_plot1 + geom_point(shape=21,color="white",aes(fill=Group),size=5) + theme_bw()
#Ellips
frog_site_clr_plot3<- frog_site_clr_plot2 + geom_path(data = x1$df_ellipse, aes(x=x, y=y, color=Group), show.legend = FALSE)
#Colour
frog_site_clr_plot4<- frog_site_clr_plot3 + labs(x="PC1 (10.97%)",y="PC2 (7.29%)",fill="Site") + plotopts + scale_fill_manual(values=sitecols) + scale_color_manual(values=sitecols) + theme(legend.position = "bottom")
frog_site_clr_plot4
ggsave2('CLR Beta Diversity By SITE.pdf',frog_site_clr_plot4,width=22,height=15,units="cm")1.5 BETA DIVERSITY
1.5.1 Barplots of Community Structure
We can summarise our microbial communities using compositional barplots, which provide a quick and accessible way of visualising differences in taxonomy at different hierarchies.
The first step here is to remove samples in our phyloseq object with low reads, otherwise this will break our averaging
#Filter Samples with no Data - FROGS
ranoides_ps_spatial<-prune_samples(sample_sums(ps_rare)>0& sample_data(ps_rare)$study=="spatial",ps_rare)1.5.2 Top N Approach
The plots can get quite messy - so it’s not uncommon to just subset to the top ‘N’ of a taxonomic group for easier visualisation
#What Are the Names of the most abundant phyla?
ranoides_ps_spatial_phylum<- ranoides_ps_spatial %>% aggregate_taxa(level="Phylum")
physeq_top5phyla = names(sort(taxa_sums(ranoides_ps_spatial_phylum), TRUE)[1:5])
physeq_top5phyla ## [1] "Proteobacteria" "Bacteroidota" "Firmicutes" "Actinobacteriota"
## [5] "Acidobacteriota"#Subset the phyloseq object to those phyla
ranoides_ps_spatial_phylumtop5<-subset_taxa(ranoides_ps_spatial,Phylum %in% physeq_top5phyla)
sample_data(ranoides_ps_spatial_phylumtop5)$site_name_short<-sapply(strsplit(sample_data(ranoides_ps_spatial_phylumtop5)$site,"-"),"[",1)
#Remake Our Graph but with no grouping (samples)
ranoides_ps_spatial_phylumplot <- ranoides_ps_spatial_phylumtop5 %>%
aggregate_taxa(level = "Phylum") %>%
microbiome::transform(transform = "compositional") %>%
plot_composition(sample.sort = "Proteobacteria",group_by="site_name_short") + theme_bw(base_size=15) + scale_fill_manual(values = moma.colors("Levine2", direction=-1)) + labs(fill="Phylum",x="") + theme(legend.position = "top",axis.text.x = element_blank(),
axis.ticks.x = element_blank())
ranoides_ps_spatial_phylumplot 
1.6 Combined Spatial Fig.
spatial_fig<-plot_grid(site_new_richness_bayes2,frog_site_clr_plot4,labels="AUTO",label_size = 25)
spatial_fig2<-plot_grid(spatial_fig,ranoides_ps_spatial_phylumplot,nrow=2,labels=c("","C"),label_size = 25)
spatial_fig2
ggsave2('Fig 3 Spatial Microbiota Patterns.pdf',width=30,height=25,units="cm")
ggsave2('Fig 3 Spatial Microbiota Patterns.png',width=30,height=25,units="cm")1.7 Heatmaps
1.7.1 Heatmaps in phyloseq
##### HEATMAP (baked into phyloseq)
#Subset to Most 50 abundant OTUs
ps_rare_top50 <- prune_taxa(names(sort(taxa_sums(ps_rare),TRUE)[1:50]), ps_rare)
ps_rare_top50_spatial<-prune_samples(sample_data(ps_rare_top50)$study=="spatial",ps_rare_top50)
#Plot Heatmap with X axis ordered Inter-Sample Distance
plot_heatmap(ps_rare_top50_spatial,"NMDS",distance = "bray") ## Warning in scale_fill_gradient(low = low, high = high, trans = trans, na.value
## = na.value): log-4 transformation introduced infinite values.
#And again, with explicit ordering of samples by SITE
plot_heatmap(ps_rare_top50_spatial,"NMDS",distance = "bray",sample.label="site") ## Warning in scale_fill_gradient(low = low, high = high, trans = trans, na.value
## = na.value): log-4 transformation introduced infinite values.
1.7.2 Other Heatmaps (pheatmap)
#Generate Metadata for Plotting Colours (SITE ID)
site_data<-data.frame(Site=sample_data(ps_rare_top50_spatial)$site)
rownames(site_data)<-rownames(sample_data(ps_rare_top50_spatial))
#Strip Out OTU_Table
ps_rare_top50_spatial_otu<-otu_table(ps_rare_top50_spatial)
#Plot Heatmap - No Explicit Clustering
pheatmap(ps_rare_top50_spatial_otu,cluster_cols = FALSE,scale="row",annotation_col = site_data,border_color = "white")
Now cluster columns by similarity
#As Above, But Order Columns (Samples) by Similarity
pheatmap(ps_rare_top50_spatial_otu,cluster_cols = TRUE,scale="row",annotation_col = site_data,border_color = "white")
1.8 PERMANOVA: Statistical Testing of Beta Diversity
1.8.1 Getting the Data Ready
#Convenience Functions
#Function to Extract the Data from phyloseq in a format vegan can understand
vegan_otu <- function(physeq) {
OTU <- otu_table(physeq)
if (taxa_are_rows(OTU)) {
OTU <- t(OTU)
}
return(as(OTU, "matrix"))
}
#Extract Matrix and Sample Data
frog_spatial_clr_v<-vegan_otu(ps_spatial_clr)
frog_spatial_clr_s<-as(sample_data(ps_spatial_clr),"data.frame")
frog_spatial_clr_v2<-frog_spatial_clr_v[-which(is.na(frog_spatial_clr_s$sv_mm)),]
frog_spatial_clr_s2<-frog_spatial_clr_s[-which(is.na(frog_spatial_clr_s$sv_mm)),]
## Perform PERMANOVA
frog_clr_perm <- adonis2(frog_spatial_clr_v2 ~ site_name_short + sv_mm ,data=frog_spatial_clr_s2, permutations=999,method="euclidean",by="terms")
frog_clr_perm## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = frog_spatial_clr_v2 ~ site_name_short + sv_mm, data = frog_spatial_clr_s2, permutations = 999, method = "euclidean", by = "terms")
## Df SumOfSqs R2 F Pr(>F)
## site_name_short 2 9403 0.07058 1.6777 0.001 ***
## sv_mm 1 3329 0.02499 1.1880 0.125
## Residual 43 120495 0.90444
## Total 46 133227 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## Perform PERMANOVA
frog_clr_perm <- adonis2(frog_spatial_clr_v2 ~ site + sv_mm ,data=frog_spatial_clr_s2, permutations=999,method="euclidean",by="mar")
frog_clr_perm ## Permutation test for adonis under reduced model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = frog_spatial_clr_v2 ~ site + sv_mm, data = frog_spatial_clr_s2, permutations = 999, method = "euclidean", by = "mar")
## Df SumOfSqs R2 F Pr(>F)
## site 5 19543 0.14669 1.4168 0.002 **
## sv_mm 1 3291 0.02470 1.1928 0.114
## Residual 40 110348 0.82827
## Total 46 133227 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 frog_clr_ph_perm <- adonis2(frog_spatial_clr_v2 ~ p_h_mean ,data=frog_spatial_clr_s2, permutations=999,method="euclidean",by="mar")
frog_clr_ph_perm ## Permutation test for adonis under reduced model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = frog_spatial_clr_v2 ~ p_h_mean, data = frog_spatial_clr_s2, permutations = 999, method = "euclidean", by = "mar")
## Df SumOfSqs R2 F Pr(>F)
## p_h_mean 1 4107 0.03082 1.4312 0.028 *
## Residual 45 129120 0.96918
## Total 46 133227 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 12 GLLVM
## OTUs Top 50
ranoides_genus<-aggregate_top_taxa2(ps_clean_spatial, "Genus", top = 50)
ranoides_clr_scaled <-microbiome::transform(ranoides_genus, transform = "clr")
sample_data(ranoides_clr_scaled)$site_name_short<-sapply(strsplit(sample_data(ranoides_clr_scaled)$site,"-"),"[",1)
#Extract
ys <- data.frame(t(otu_table(ranoides_clr_scaled)))
names(ys) <-taxa_names(ranoides_clr_scaled)
#Predictors
Xs<-data.frame(sample_data(ranoides_clr_scaled)) %>% dplyr::select(site_name_short,p_h_mean,sv_mm)
#Strip Missing SVL
sv_miss<-which(is.na(Xs$sv_mm))
Xs<-Xs[-sv_miss,]
ys<-ys[-sv_miss,]
## Model
# fit_reduced_scaled <- gllvm(ys, Xs,
# num.lv = 2,
# formula = ~ p_h_mean + sv_mm,
# family = "gaussian",starting.val='zero',
# row.eff = "random",randomX=~site_name_short)
#
fit_reduced_scaled <- gllvm(ys, Xs,
num.lv = 2,
formula = ~ site_name_short + sv_mm,
family = "gaussian",starting.val='zero')## Warning in gllvm(ys, Xs, num.lv = 2, formula = ~site_name_short + sv_mm, : There are rows full of zeros in y. coefplot(fit_reduced_scaled)
#plot(fit_reduced_scaled)
#Estimates
df<-coef(fit_reduced_scaled)
est_df<-data.frame(df$Intercept)
est_df2<-data.frame(df$Xcoef)
est_df3<-merge(est_df, est_df2, by = 0)
#Order genera
row.names(est_df3)<-est_df3$Row.names
est_df3<-est_df3[colnames(ys),]
names(est_df3)[1]<- "Genus"
names(est_df3)[2]<- "Intercept"
names(est_df3)[3:5]<-c("Site M", "Site P", "SVL")
### COnfidence Intervals
confint_df<-data.frame(confint(fit_reduced_scaled))
###Strip Out Individual Datasets
##Site M
siteM_locs<-grepl("Xcoef.site_name_shortM",rownames(confint_df))
siteM_data<-cbind(est_df3[,c("Genus","Site M")],confint_df[siteM_locs,])
colnames(siteM_data)<-c("Genus","Estimate","l95","u95")
siteM_data$trait<-"Site M"
## Site P
siteP_locs<-grepl("Xcoef.site_name_shortP",rownames(confint_df))
siteP_data<-cbind(est_df3[,c("Genus","Site P")],confint_df[siteP_locs,])
colnames(siteP_data)<-c("Genus","Estimate","l95","u95")
siteP_data$trait<-"Site P"
##SVL
svl_locs<-grepl("Xcoef.sv_mm",rownames(confint_df))
svl_data<-cbind(est_df3[,c("Genus","SVL")],confint_df[svl_locs,])
colnames(svl_data)<-c("Genus","Estimate","l95","u95")
svl_data$trait<-"SVL"
ranoides_mod_plotdata<-rbind(siteM_data,siteP_data,svl_data)
#Strip Out Overly Long Multi-Genus Assignment
# ranoides_mod_plotdata2<-ranoides_mod_plotdata[-which(grepl("Allorhizobium-",ranoides_mod_plotdata$Genus)),]
ranoides_mod_plotdata3<- ranoides_mod_plotdata %>% group_by(trait) %>% arrange(Estimate,.by_group=T)
ranoides_mod_plotdata3$Genus<-factor(ranoides_mod_plotdata3$Genus,levels=unique(ranoides_mod_plotdata3$Genus))
#Significance
ranoides_mod_plotdata3$Sig<- !data.table::between(0, ranoides_mod_plotdata3$l95, ranoides_mod_plotdata3$u95)
#Filter by Significance
ranoides_mod_plotdata_sig<-subset(ranoides_mod_plotdata3,Sig==T)
### plot
ranoides_gllvm_plot1<-ggplot(ranoides_mod_plotdata_sig,aes(x=Estimate,y=Genus)) + geom_errorbarh(aes(xmin=l95,xmax=u95),color="grey22") + geom_point(size=5,shape=21,color="gray40",fill="cornflowerblue",alpha=1) + theme_bw()
ranoides_gllvm_plot2<- ranoides_gllvm_plot1 + theme_bw(base_size = 15) + geom_vline(xintercept=0,linetype="dashed") + guides(fill="none",color="none") + facet_wrap(.~trait,scales = "free_x")
ranoides_gllvm_plot2
ggsave2('Ranoides GLLVM Output.png',ranoides_gllvm_plot2,width=20,height=15,units="cm")3 UpSet Plot
#Create A Presence/Absence Dataframe for UpSet Plots
upsetda <- MicrobiotaProcess::get_upset(obj=ps_clean_spatial, factorNames="site_name_short")
#Copy Across Phylum Level Data
Phylumvals<-as.vector(tax_table(ps_clean_spatial)[,2][match(rownames(upsetda),rownames(tax_table(ps_clean_spatial)))])
#Tabluate and Rank by Frequency of Phylum
top5phyla<- data.frame(table(Phylumvals)) %>% arrange(desc(Freq))
#Overwrite any not in the top X as "Other"
Phylumvals[which(!Phylumvals %in% as.character(top5phyla$Phylumvals[1:5]))]<-"Other"
#Copy Across
upsetda$Phylum<-factor(Phylumvals,levels=c("Actinobacteriota","Bacteroidota","Firmicutes","Proteobacteria","Verrucomicrobiota","Other"))
#install.packages("ComplexUpset")
#library(ComplexUpset)
#Sets
sites<-colnames(upsetda)[1:3]
#Trim ASVs not in Any of the 6
upsetda_trimmed<-upsetda[apply(upsetda[,1:3],1,sum)>0,]
#Basic Upset Plot
upset(upsetda_trimmed,sites)
#Cols to Overlap With Phylum Plot
upsetcols<-c(moma.colors("Levine2",direction = -1)[c(2:5)],moma.colors("Rattner",3)[2:3])
#Add a Stacked Barplot for Phylum
ranoides_upset<-upset(
upsetda_trimmed,sites,set_size=FALSE,
themes=upset_default_themes(text=element_text(size=16),legend.position="top"),base_annotations=list(
'Intersection size'=intersection_size(counts=TRUE)),
annotations = list(
'Phylum'=(
ggplot(mapping=aes(fill=Phylum))
+ geom_bar(stat='count', position='fill')
+ scale_y_continuous(labels=scales::percent_format())
+ ylab('Phylum') + scale_fill_manual(values = upsetcols)
)
),
width_ratio=0.1
)
ranoides_upset
cowplot::ggsave2('Ranoides Upset Plot and PhylumStack.pdf',ranoides_upset,width=20,height=25,units="cm")
cowplot::ggsave2('Ranoides Upset Plot and PhylumStack.png',ranoides_upset,width=20,height=25,units="cm")4 Sample Similarity Network
library(NetCoMi)## Loading required package: SpiecEasi##
## Attaching package: 'SpiecEasi'## The following object is masked from 'package:MASS':
##
## fitdistr## library(SpiecEasi)
#Build Network
spatial_sample_net <- netConstruct(ps_spatial_clr,,dataType = "counts",measure = "aitchison",
measurePar = list(method = "mb",
pulsar.params = list(rep.num = 10),
symBetaMode = "ave"),
normMethod = "clr",
sparsMethod = "knn")## Checking input arguments ...## Done.## Infos about changed arguments:## Zero replacement needed for measure "aitchison". "multRepl" used.
## Counts normalized to fractions for measure "aitchison".## 1038 taxa and 48 samples remaining.##
## Zero treatment:## Data contains no zeros.##
## Normalization:## Counts normalized by total sum scaling.##
## Calculate 'aitchison' dissimilarities ... Done.
##
## Sparsify dissimilarities via 'knn' ... Done.#Analyse and Plot
spatial_sample_net_an <- netAnalyze(spatial_sample_net, clustMethod = "cluster_fast_greedy")## Warning: The `scale` argument of `eigen_centrality()` always as if TRUE as of igraph
## 2.1.1.
## ℹ Normalization is always performed
## ℹ The deprecated feature was likely used in the NetCoMi package.
## Please report the issue at <https://github.com/stefpeschel/NetCoMi/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
######### NAIVE PLOT WITH STATISTICAL CLUSTERS
plot(spatial_sample_net_an,rmSingles = T)
######## VERSION TO MATCH ASV NETWORK - CLUSTERS ARE COLOURED, SAMPLE TYPE = SHAPE
#Features to Show Sample Type (ALSO USED BELOW)
featvec<-frog_spatial_clr_s$site_name_short
names(featvec)<-frog_spatial_clr_s$sample
#Site Labels
site_labs<-frog_spatial_clr_s$site_name_short
names(site_labs)<-frog_spatial_clr_s$sample
png('Ranoides Sample Network ClusterCol.png',res=600,width=40,height=30,units="cm")
plot(spatial_sample_net_an,rmSingles = T,
repulsion=1.2,
labelScale = FALSE,
cexLabels = 1,
cexNodes = 1,layout="spring",nodeSize = "fix",
featVecShape=featvec,nodeColor = "cluster",featVecCol =featvec,labels=site_labs,colorVec = brewer.pal(6,"Dark2")[1:3])
legend("topright", title = "",legend=c("BdPG","M","P"), bty = "n",pch=c(16,15,17),xpd=T,cex = 3, pt.cex = 5,y.intersp=1,x.intersp = 1,col="gray")
dev.off() ## quartz_off_screen
## 25 TEMPORAL ANALYSIS
5.1 Translocation Frog Mass Dynamics
transmass<-read.csv('Ranoides Trip 2 data.csv')
transmass<-janitor::clean_names(transmass)
#Pivot to Long Format
library(tidyr)
transmass$weight_post[transmass$weight_post=="No recapture"]<-NA
transmass$weight_post<-as.numeric(transmass$weight_post)
tmass_long<-pivot_longer(transmass,cols = starts_with("weight"),names_to="time")
tmass_long$time<-factor(tmass_long$time,levels=c("weight_pre","weight_post"),labels = c("PreTranslocation","PostTranslocation"))
tmass_long$frog_number<-factor(tmass_long$frog_number)5.2 Richness over Time
frog_temporal<-subset(frog_rich,study=="temporal")
frog_temporal$timepoint<-factor(frog_temporal$timepoint,levels=c("PreTranslocation","PostTranslocation"))
frog_temporal<-frog_temporal %>% arrange(id)
frog_temporal$id<-factor(frog_temporal$id)
frog_temporal$sex<-transmass$sex[match(frog_temporal$id,transmass$frog_number)]
#Plotting Colours
library(Polychrome)
IDcols<-data.frame(id=seq(14),colval = palette36.colors(14))
frog_temporal$colval<-IDcols$colval[match(frog_temporal$id,IDcols$id)]
##Plot
frog_temporal_richplot1<-ggplot(frog_temporal,aes(x=timepoint,y=Observed)) + geom_line(aes(group=id,color=id),position = position_dodge(width=0.2),linewidth=1.2) + geom_point(size=6,aes(fill=id,shape=sex),colour="white",position = position_dodge(width=0.2)) + theme_bw()
frog_temporal_richplot2<- frog_temporal_richplot1 + plotopts + labs(x="Time Point",y="Alpha Diversity \n(Richness)") + guides(fill="none",color="none") + scale_fill_manual(values=IDcols$colval) + scale_color_manual(values=IDcols$colval) + scale_shape_manual(values=c(21,23)) + guides(shape = guide_legend(override.aes = list(size = 5,shape=c(21,23),fill="grey")),fill=guide_legend(override.aes = list(shape=21))) + labs(shape="Sex",fill="Frog ID") + theme(axis.text.x=element_text(size=12))
frog_temporal_richplot2
# frog_temporal_richplot3<-ggplot(frog_temporal,aes(x=timepoint,y=Observed)) + geom_line(aes(group=id,color=id),position = position_dodge(width=0.2),linewidth=1.2) + geom_point(size=6,shape=21,aes(fill=id),colour="white",position = position_dodge(width=0.2)) + theme_bw()
# frog_temporal_richplot4<- frog_temporal_richplot1 + plotopts + labs(x="Time Point",y="Alpha Diversity \n(Richness)") + guides(fill="none",color="none") + scale_fill_manual(values=unname(palette36.colors(14))) + scale_color_manual(values=unname(palette36.colors(14)))
# frog_temporal_richplot4#### PLOTS
transloccols<-IDcols$colval[match(tmass_long$frog_number,IDcols$id)]
# ggplot(tmass_long,aes(x=time,y=value)) + geom_point(aes(fill=frog_number),shape=21,size=5) + theme_bw(base_size=15) + labs(y="Mass (g)",x="Time Point") + guides(fill="none") + scale_fill_manual(values=IDcols) + scale_color_manual(values=transloccols)
#
massplot1<-ggplot(tmass_long,aes(x=time,y=value)) + geom_line(aes(group=frog_number,color=frog_number),position=position_dodge(width=0.2)) + geom_point(aes(fill=frog_number,shape=sex),position=position_dodge(width=0.2),size=7,color="white") + theme_bw(base_size=20) + labs(y="Mass (g)",x="Time Point") + guides(color="none") + scale_fill_manual(values=unname(palette36.colors(14))) + scale_color_manual(values=unname(palette36.colors(14))) + scale_shape_manual(values=c(21,23))
massplot2<- massplot1 + guides(shape = guide_legend(override.aes = list(size = 5,shape=c(21,23),fill="grey")),fill=guide_legend(override.aes = list(shape=21))) + labs(shape="Sex",fill="Frog ID")
massplot2## Warning: Removed 5 rows containing missing values or values outside the scale range
## (`geom_line()`).## Warning: Removed 5 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave2('Translocation Mass Dynamics.png',massplot2,width=10,height=10)## Warning: Removed 5 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 5 rows containing missing values or values outside the scale range
## (`geom_point()`).5.3 Plot Beta Diversity Over Time
5.4 Beta Diversity
### Just Frogs
ps_clean_temporal<-prune_samples(sample_data(ps_clean_samples)$study=="temporal",ps_clean_samples)
#CLR TRANSFORM
ps_temporal_clr<-microbiome::transform(ps_clean_temporal, 'clr')
#Ordinate
ord_temporal_clr <- phyloseq::ordinate(ps_temporal_clr, "RDA")
#phyloseq::plot_ordination(ps_spatial_clr, ord_spatial_clr, type="samples", color="site_name")
#Extract Scores
temporal_scores<-as.data.frame(scores(ord_temporal_clr)$sites)
temporal_scores$sample<-rownames(temporal_scores)
temporal_scores<-left_join(temporal_scores,frog_meta,"sample")
temporal_scores$id<-factor(temporal_scores$id)
temporal_scores$timepoint<-factor(temporal_scores$timepoint,levels=c("PreTranslocation","PostTranslocation"))5.4.1 Ordination Over TIme
#GG Version
gg_ordiplot(ord_temporal_clr,sample_data(ps_temporal_clr)$timepoint,spiders = TRUE,ellipse = T)
#### Extract for Plotting
#Plot and Extract for GGPLot2
time1 <- gg_ordiplot(ord_temporal_clr,sample_data(ps_temporal_clr)$timepoint,spiders = TRUE,ellipse = T,plot=F)
time1_plotdata<- time1$df_ord
time1_plotdata$Group<-factor(time1_plotdata$Group,levels=c("PreTranslocation","PostTranslocation"))
time1$df_spiders$Group<-factor(time1$df_spiders$Group,levels=c("PreTranslocation","PostTranslocation"))
#Plot Spiders
frog_time_clr_plot1<-ggplot(time1_plotdata,aes(x=x,y=y)) + geom_segment(data=time1$df_spiders, aes(x=cntr.x, xend=x, y=cntr.y, yend=y, color=Group), show.legend = FALSE)
#Plot Points
frog_time_clr_plot2<- frog_time_clr_plot1 + geom_point(shape=21,color="white",aes(fill=Group),size=5) + theme_bw()
#Ellips
frog_time_clr_plot3<- frog_time_clr_plot2 + geom_path(data = time1$df_ellipse, aes(x=x, y=y, color=Group), show.legend = FALSE)
#Colour
frog_time_clr_plot4<- frog_time_clr_plot3 + labs(x="PC1 (35.27%)",y="PC2 (7.83%)",fill="Time") + plotopts + scale_fill_brewer(palette = "Set2") + scale_color_brewer(palette = "Set2") + theme(legend.position = "top") + guides(fill=guide_legend(nrow=2,byrow=TRUE))
frog_time_clr_plot4
ggsave2('CLR Beta Diversity Translocation.pdf',frog_time_clr_plot4,width=18,height=15,units="cm")
###### PERMANOVA
#Extract Matrix and Sample Data
frog_temporal_clr_v<-vegan_otu(ps_temporal_clr)
frog_temporal_clr_s<-as(sample_data(ps_temporal_clr),"data.frame")
frog_temporal_clr_s$id<-factor(frog_temporal_clr_s$id)
frog_temporal_clr_s$id_num<-as.numeric(frog_temporal_clr_s$id)
#Permutation Structure
library("permute")
h <- how(plots = Plots(strata = frog_temporal_clr_s$id, type = "none"),
nperm = 999,within = Within(type="series"))
# Run PERMANOVA
adonis2(frog_temporal_clr_v ~ timepoint + id, data=frog_temporal_clr_s,method="euclidean",permutations = h,by="mar")## 'nperm' >= set of all permutations: complete enumeration.## Set of permutations < 'minperm'. Generating entire set.## Permutation test for adonis under reduced model
## Marginal effects of terms
## Plots: frog_temporal_clr_s$id, plot permutation: none
## Permutation: series
## Number of permutations: 255
##
## adonis2(formula = frog_temporal_clr_v ~ timepoint + id, data = frog_temporal_clr_s, permutations = h, method = "euclidean", by = "mar")
## Df SumOfSqs R2 F Pr(>F)
## timepoint 1 2556 0.03248 0.5104 0.7422
## id 13 38638 0.49103 0.5935 0.7031
## Residual 7 35053 0.44547
## Total 21 78688 1.000005.4.2 Procrustes Rotation
#Make Sure Only Samples is Measured Twice Occur
temporal_idtab<-table(sample_data(ps_clean_temporal)$id)
temporal_idtab_rep2<-names(temporal_idtab)[temporal_idtab>1]
#Time Point 1
temporal_pretrans<-prune_samples(sample_data(ps_clean_temporal)$id %in% temporal_idtab_rep2 & sample_data(ps_clean_temporal)$timepoint=="PreTranslocation",ps_clean_temporal)
temporal_pretrans_clr<-microbiome::transform(temporal_pretrans, 'clr')
ord_temporal_pretrans_clr <- phyloseq::ordinate(temporal_pretrans_clr, "RDA")
#Time Point 2
temporal_posttrans<-prune_samples(sample_data(ps_clean_temporal)$id %in% temporal_idtab_rep2 & sample_data(ps_clean_temporal)$timepoint=="PostTranslocation",ps_clean_temporal)
temporal_posttrans_clr<-microbiome::transform(temporal_posttrans, 'clr')
ord_temporal_posttrans_clr <- phyloseq::ordinate(temporal_posttrans_clr, "RDA")
## RUN PROCRUSTES
translocation_procrustes<-procrustes(X=ord_temporal_posttrans_clr,Y=ord_temporal_pretrans_clr,symmetric = T)
plot(translocation_procrustes)
protest(X=ord_temporal_posttrans_clr,Y=ord_temporal_pretrans_clr,scores="sites")##
## Call:
## protest(X = ord_temporal_posttrans_clr, Y = ord_temporal_pretrans_clr, scores = "sites")
##
## Procrustes Sum of Squares (m12 squared): 0.7602
## Correlation in a symmetric Procrustes rotation: 0.4897
## Significance: 0.381
##
## Permutation: free
## Number of permutations: 9995.4.3 Temporal Network
library(NetCoMi)
library(SpiecEasi)
#Build Network
temporal_sample_net <- netConstruct(ps_temporal_clr,dataType = "counts",measure = "aitchison",
measurePar = list(method = "mb",
pulsar.params = list(rep.num = 10),
symBetaMode = "ave"),
normMethod = "clr",
sparsMethod = "knn")## Checking input arguments ... Done.
## Infos about changed arguments:
## Zero replacement needed for measure "aitchison". "multRepl" used.
## Counts normalized to fractions for measure "aitchison".
##
## 1038 taxa and 22 samples remaining.
##
## Zero treatment:
## Data contains no zeros.
##
## Normalization:
## Counts normalized by total sum scaling.
##
## Calculate 'aitchison' dissimilarities ... Done.
##
## Sparsify dissimilarities via 'knn' ... Done.#Analyse and Plot
temporal_sample_net_an <- netAnalyze(temporal_sample_net, clustMethod = "cluster_fast_greedy")
######### NAIVE PLOT WITH STATISTICAL CLUSTERS
plot(temporal_sample_net_an,rmSingles = T)
######## VERSION TO MATCH ASV NETWORK - CLUSTERS ARE COLOURED, SAMPLE TYPE = SHAPE
#Features to Show Sample Type (ALSO USED BELOW)
featvec_temp<-frog_temporal_clr_s$timepoint
names(featvec_temp)<-frog_temporal_clr_s$sample
#Site Labels
id_labs<-as.character(frog_temporal_clr_s$id)
names(id_labs)<-frog_temporal_clr_s$sample
png('Ranoides Temporal Network ClusterCol.png',res=600,width=40,height=30,units="cm")
plot(temporal_sample_net_an,rmSingles = T,
repulsion=1.2,
labelScale = FALSE,
cexLabels = 3,
cexNodes = 1,layout="spring",nodeSize = "fix",
featVecShape=featvec_temp,nodeColor = "cluster",labels=id_labs,colorVec = brewer.pal(6,"Dark2")[1:3])
legend("topright", title = "",legend=c("PreTranslocation","PostTranslocation"), bty = "n",pch=c(15,16),xpd=T,cex = 3, pt.cex = 5,y.intersp=1,x.intersp = 1,col="gray")
dev.off() ## quartz_off_screen
## 2###Heatmap
#Subset to Most 50 abundant OTUs
ps_temporal_top50 <- prune_taxa(names(sort(taxa_sums(ps_clean_temporal),TRUE)[1:50]), ps_clean_temporal)
## Only Repeat Sampled Individuals
repeat_frogs<-table(sample_data(ps_temporal_top50)$id)
ps_temporal_top50_repeatsampled <- prune_samples(sample_data(ps_temporal_top50)$id %in% names(repeat_frogs)[repeat_frogs==2], ps_temporal_top50)
#Generate Metadata for Plotting Colours (SITE ID)
id_data<-data.frame(ID=factor(sample_data(ps_temporal_top50_repeatsampled)$id),Timepoint=sample_data(ps_temporal_top50_repeatsampled)$timepoint)
rownames(id_data)<-rownames(sample_data(ps_temporal_top50_repeatsampled))
#Strip Out OTU_Table
ps_rare_top50_temporal_otu<-otu_table(rarefy_even_depth(ps_temporal_top50_repeatsampled))## You set `rngseed` to FALSE. Make sure you've set & recorded
## the random seed of your session for reproducibility.
## See `?set.seed`## ... #Order By Sample ID
sampids<-sample_data(ps_temporal_top50_repeatsampled)$sample
sampids_sequence<-rep(1,length(sampids))
sampids_sequence[grep("_2",sampids)]<-2
ps_rare_top50_temporal_otu_ordered<-ps_rare_top50_temporal_otu[, order(sample_data(ps_temporal_top50_repeatsampled)$id,sampids_sequence)]
#Set2cols
set2<-brewer.pal(5,"Set2")
ann_colors = list(
Timepoint = c(PreTranslocation=set2[1],PostTranslocation=set2[2]),
ID = palette36.colors(14))
names(ann_colors[[2]])<-seq(14)
#Rownames
rowlabels_genus<-data.frame(Genus=tax_table(ps_temporal_top50_repeatsampled)[,6],Family=tax_table(ps_temporal_top50_repeatsampled)[,5])
gen_na<-which(is.na(rowlabels_genus$Genus))
rowlabels_genus$Tax_combined<-rowlabels_genus$Genus
rowlabels_genus$Tax_combined[gen_na]<-paste0("[",rowlabels_genus$Family[gen_na],"]")
rownames(rowlabels_genus)<-rownames(ps_rare_top50_temporal_otu_ordered)
#Plot Heatmap - With Hierarchical Clustering
heatmap1<-pheatmap(ps_rare_top50_temporal_otu_ordered,cluster_cols = TRUE,scale="row",annotation_col = id_data,border_color = "white",annotation_colors = ann_colors,labels_row = rowlabels_genus[,3],show_colnames=F)
5.5 Grid Plot Individual Level Temporal
temp_grid1<-plot_grid(frog_temporal_richplot2,frog_time_clr_plot4,labels="AUTO",label_size = 20)
temp_grid2<-plot_grid(temp_grid1,heatmap1[[4]],nrow=2,labels=c("","C"),label_size = 20)
temp_grid2
ggsave2('Individual Temporal Dynamics Grid.pdf',temp_grid2,width=30,height=30,units="cm")
ggsave2('Individual Temporal Dynamics Grid.png',temp_grid2,width=30,height=30,units="cm")6 Temporal GLLVM
## OTUs Top 50
ranoides_temporal_genus<-aggregate_top_taxa2(ps_clean_temporal, "Genus", top = 50)
ranoides_temporal_clr_scaled <-microbiome::transform(ranoides_temporal_genus, transform = "clr")
#Extract
temporal_ys <- data.frame(t(otu_table(ranoides_temporal_clr_scaled)))
names(temporal_ys) <-taxa_names(ranoides_temporal_clr_scaled)
#Predictors
temporal_Xs<-data.frame(sample_data(ranoides_temporal_clr_scaled)) %>% dplyr::select(timepoint,id)
temporal_Xs$timepoint<-factor(temporal_Xs$timepoint,levels=c("PreTranslocation","PostTranslocation"))
## Model
fit_temporal1 <- gllvm(temporal_ys, temporal_Xs,
num.lv = 2,
formula = ~ timepoint,
family = "gaussian",starting.val='zero',
row.eff = "random",randomX=~id)
coefplot(fit_temporal1)
#Estimates
dftemp<-coef(fit_temporal1)
est_dftemp<-data.frame(dftemp$Intercept)
est_dftemp2<-data.frame(dftemp$Xcoef)
est_dftemp3<-merge(est_dftemp, est_dftemp2, by = 0)
#Order genera
row.names(est_dftemp3)<-est_dftemp3$Row.names
est_dftemp3<-est_dftemp3[colnames(temporal_ys),]
names(est_dftemp3)[1]<- "Genus"
names(est_dftemp3)[2]<- "Intercept"
### COnfidence Intervals
confint_dftemp<-data.frame(confint(fit_temporal1))
###Strip Out Individual Datasets
##TimePoint
posttrans_main<-grepl("^Xcoef.timepointPostTranslocation",rownames(confint_dftemp))
post_maindata<-cbind(est_dftemp3[,c("Genus","timepointPostTranslocation")],confint_dftemp[posttrans_main,])
colnames(post_maindata)<-c("Genus","Estimate","l95","u95")
## Factor
post_maindata2<- post_maindata %>% arrange(Estimate,.by_group=T)
post_maindata2$Genus<-factor(post_maindata2$Genus,levels=unique(post_maindata2$Genus))
#Significance
post_maindata2$Sig<- !data.table::between(0, post_maindata2$l95, post_maindata2$u95)
sig_col<-moma.colors("Andri",3)[3]
#Filter by Sig
post_maindata2_sig<-subset(post_maindata2,Sig==T)
### plot
ranoides_temporal_gllvm_plot1<-ggplot(post_maindata2_sig,aes(x=Estimate,y=Genus)) + geom_errorbarh(aes(xmin=l95,xmax=u95),color="gray22",alpha=1,height=0.15) + geom_point(size=5,shape=21,fill="cornflowerblue",alpha=1) + theme_bw()
ranoides_temporal_gllvm_plot2<- ranoides_temporal_gllvm_plot1 + theme_bw(base_size = 15) + geom_vline(xintercept=0,linetype="dashed") + scale_color_manual(values=c("gray40",sig_col)) + scale_fill_manual(values=c("white",sig_col)) + guides(fill="none",color="none")
ranoides_temporal_gllvm_plot2
ggsave2('Ranoides Temporal GLLVM Output.png',ranoides_temporal_gllvm_plot2,width=15,height=15,units="cm")
#### Correlations
cr1<-data.frame(getResidualCor(fit_reduced_scaled))#
names(cr1)<-names(ys)
names(cr1)<-abbreviate(names(cr1), minlength = 15)
rownames(cr1)<-abbreviate(rownames(cr1), minlength = 15 )
library(rstatix)##
## Attaching package: 'rstatix'## The following object is masked from 'package:janitor':
##
## make_clean_names## The following object is masked from 'package:MASS':
##
## select## The following object is masked from 'package:stats':
##
## filter #devtools::install_github("kassambara/ggcorrplot")
library(ggcorrplot)##
## Attaching package: 'ggcorrplot'## The following object is masked from 'package:rstatix':
##
## cor_pmat #install.packages("ggpubr")
library(ggpubr)##
## Attaching package: 'ggpubr'## The following object is masked from 'package:cowplot':
##
## get_legend cr2<-cor_pmat(cr1)7 FUNCTIONAL ANALYSIS
7.1 TEMPORAL Prepare for MicFunPred
Adapted from https://forum.qiime2.org/t/importing-dada2-and-phyloseq-objects-to-qiime-2/4683
## Trim Newt1
# ranoides_filter<-prune_taxa(taxa_sums(ps_clean_temporal)>200,ps_clean_temporal)
#
# #Make Fasta File From Sequences
# #remotes::install_github("alexpiper/seqateurs")
# library(seqateurs)
# ps_to_fasta(ranoides_filter)
#
# #Export Feature Table
# #library(biomformat);packageVersion("biomformat")
# otu<-(as(otu_table(ranoides_filter),"matrix")) # 't' to transform if taxa_are_rows=FALSE
# #Rename ASVs to Match FASTA above
# rownames(otu)<-paste0("ASV_",seq(nrow(otu)))
# #Fix Sample Names
# true_samples_ranoides_filter<-colnames(otu)
# colnames(otu)<-paste0("sample",seq(ncol(otu)))
# #Write File
# # otu_biom<-make_biom(data=otu)
# # write_biom(otu_biom,"Newt1_filter.biom")
# # write.table(otu,'ranoides_temporal_filter.tsv',sep="\t",row.names=TRUE, col.names=TRUE, quote=FALSE)
## Sample Data
#ranoides_filter_samplemeta<-as(sample_data(ranoides_filter),"data.frame")#
#write.csv(ranoides_filter_samplemeta,"Ranoides_Temporal_Filter_sample-metadata.csv")7.2 Predicted Function Analysis
#KEGG pathway abundances
predfunc<-read.table('Ranoides Functional Analysis/KEGG_pathways_MinPath_prunned.tsv',sep="\t",header=T)
#Sample Data
predfunc_samples<-read.csv('Ranoides Functional Analysis/Ranoides_Temporal_Filter_sample-metadata.csv')
#Ordinate
library(vegan)
predfunc_ord<-vegan::rda(t(predfunc[,-1]))
gg_ordiplot(predfunc_ord,predfunc_samples$timepoint)
#### Extract for Plotting
#Plot and Extract for GGPLot2
timefunc1 <- gg_ordiplot(predfunc_ord,predfunc_samples$timepoint,spiders = TRUE,ellipse = T,plot=F)
timefunc1_plotdata<- timefunc1$df_ord
timefunc1_plotdata$Group<-factor(timefunc1_plotdata$Group,levels=c("PreTranslocation","PostTranslocation"))
timefunc1$df_spiders$Group<-factor(timefunc1$df_spiders$Group,levels=c("PreTranslocation","PostTranslocation"))
#Plot Spiders
frog_timefunc_clr_plot1<-ggplot(timefunc1_plotdata,aes(x=x,y=y)) + geom_segment(data=timefunc1$df_spiders, aes(x=cntr.x, xend=x, y=cntr.y, yend=y, color=Group), show.legend = FALSE)
#Plot Points
frog_timefunc_clr_plot2<- frog_timefunc_clr_plot1 + geom_point(shape=21,color="white",aes(fill=Group),size=5) + theme_bw()
#Ellips
frog_timefunc_clr_plot3<- frog_timefunc_clr_plot2 + geom_path(data = timefunc1$df_ellipse, aes(x=x, y=y, color=Group), show.legend = FALSE)
#Colour
frog_timefunc_clr_plot4<- frog_timefunc_clr_plot3 + labs(x="NMDS1",y="NMDS2",fill="Time") + plotopts + scale_fill_brewer(palette = "Set2") + scale_color_brewer(palette = "Set2") + theme(legend.position = "top")
frog_timefunc_clr_plot4
ggsave2('Functional NMDS Translocation.pdf',frog_timefunc_clr_plot4,width=18,height=15,units="cm")7.3 Functional PERMANOVA
#CLR Trasnform
temporal_func_clr <-compositions::clr(t(predfunc[,-1]))
predfunc_samples$id<-factor(predfunc_samples$id)
#Permutation Structure
library("permute")
h_func <- how(plots = Plots(strata = predfunc_samples$id, type = "none"),
nperm = 999,within = Within(type="series"))
# Run PERMANOVA
adonis2(temporal_func_clr ~ timepoint + id, data=predfunc_samples,method="euclidean",permutations = h_func,by="mar")## 'nperm' >= set of all permutations: complete enumeration.## Set of permutations < 'minperm'. Generating entire set.## Permutation test for adonis under reduced model
## Marginal effects of terms
## Plots: predfunc_samples$id, plot permutation: none
## Permutation: series
## Number of permutations: 255
##
## adonis2(formula = temporal_func_clr ~ timepoint + id, data = predfunc_samples, permutations = h_func, method = "euclidean", by = "mar")
## Df SumOfSqs R2 F Pr(>F)
## timepoint 1 258.7 0.04794 1.0047 0.4531
## id 13 3267.5 0.60561 0.9763 0.6172
## Residual 7 1802.2 0.33403
## Total 21 5395.4 1.00000